Phase diagram of silica from computer simulation 
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We evaluate the phase diagram of the "BKS" potential [Van Beest, Kramer and van Santen, Phys. 
Rev. Lett. 64, 1955 (1990)], a model of silica widely used in molecular dynamics (MD) simulations. 
We conduct MD simulations of the liquid, and three crystals (/3-quartz, coesite and stishovite) over 
wide ranges of temperature and density, and evaluate the total Gibbs free energy of each phase. The 
phase boundaries are determined by the intersection of these free energy surfaces. Not unexpectedly 
for a classical pair potential, our results reveal quantitative discrepancies between the locations of 
the BKS and real silica phase boundaries. At the same time, we find that the topology of the real 
phase diagram is reproduced, confirming that the BKS model provides a satisfactory qualitative 
description of a silica-like material. We also compare the phase boundaries with the locations of 
liquid-state thermodynamic anomalies identified in previous studies of the BKS model. 
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I. INTRODUCTION 

The melts of silica, water and a number of other molec- 
ular substances (e.g Si, Ge02, BeF2) at ambient pressure 
P form so-called "tetrahedral liquids," that is, liquids 
with properties that are strongly influenced by the oc- 
currence of a network of tetrahedrally arranged atoms. 
Such liquids display a rich spectrum of behavior, includ- 
ing density maxima [l| , dynamical anomalies 2\ , as well 
as the ability to form numerous crystal polymorphs [2j. 
Evidence also exists that liquid-liquidphase transitions 
occur in some of these systems H, M> |^ . Yet a detailed 
understanding of the commonalities among these mate- 
rials is hampered by our incomplete knowledge of their 
properties under comparable conditions. For example, 
we have extensive knowledge of liquid water for temper- 
atures T near and above the melting temperature T,„ for 
T > 0.85T,„, but the behavior below this range remains 
a subject of debate [3- Conversely, we have detailed 
knowledge of molten silica at ambient P for T < Tm , but 
a much less complete picture of the behavior at higher T 
andP 0. 

Computer simulations have contributed to filling this 
knowledge gap by providing numerical estimates of liq- 
uid behavior outside the range of current experiments. 
However, a key element has been missing from the de- 
scription of many of these model materials: their phase 
diagrams. In an experimental study of a molecular liq- 
uid, knowledge of the phase diagram - that is, the co- 
existence boundaries demarcating the stability fields of 
the liquid, gas and the various crystal phases - provides 
a vital reference that elucidates the observed thermody- 
namic, dynamic, and structural properties of the liquid 
phase. Simulations of molecular liquids are commonly 
based on semi-empirical classical interaction potentials 
that cannot be expected to precisely reproduce the ex- 



perimentally known phase diagrams of the real material. 
It is perhaps for this reason that comprehensive phase 
diagrams have not yet been developed for the simulation 
models used widely to study the complexities of impor- 
tant molecular liquids, such as water and silica. How- 
ever, as a consequence, it has not been possible to self- 
consistently relate the behavior found in simulations to 
the relevant phase boundaries of the model system, as 
would normally occur in an experimental study. 

With these motivations, we here focus on the BKS 
model of silica 9J . The BKS model has played an impor- 
tant role over the last decade in numerous studies of silica 
and related materials. For example, the BKS model has 
been used in studies of pressure induced amorphization 
of quartz [l^ , the a to (3 quartz phase transition djj, [l^ , 
the fragile-to-strong dynamical crossover in liquid sil- 
ica [i^ [ij, [13, the possibility of liquid-hquid phase 
separation in silica p\, and in the study of the generic 
topological and entropic properties of random tetrahe- 
dral networks 16]. Despite this interest in the BKS 
model, only fragments of specific crystal-crystal phase 
boundaries have been located, such as the a-[3 quartz 
transition. To our knowledge no data currently exists for 
the melting lines, though the liquid-gas coexistence curve 
has been located for a model similar to BKS jl3] • In this 
paper we report the phase diagram of the BKS model, 
finding the stability fields in the P-T plane for the liquid 
phase, and three of the prominent crystal phases of real 
silica, stishovite, coesite and /?-quartz. 



II. METHODS 

We use the BKS potential, modified at short range to 
prevent unphysical "fusion" events, and at long range to 
reduce the system size dependence of measured proper- 
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FIG. 1: Location of points in the V-T plane at which we con- 
duct simulations of the liquid (circles), stishovite (squares), 
coesite (triangles) and /3-quartz (crosses). The large symbols 
locate reference states {Vr, Tr) at which the entropy of each 
phase is evaluated directly. 
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FIG. 2: Examples of fitted and interpolated data for the liq- 
uid phase: (a) values of E along the V = 6.655 cm^ mol~^ 
isochore, fitted with a cubic polynomial (line); (b) values of 
P for V — 6.655 cm^ mol~^, fitted with a quartic polynomial 
(line); (c) interpolated values of E along the T = 4000 K 
isotherm, fitted with a cubic spline (line); (d) interpolated 
values of P for T — 4000 K, fitted with a cubic spline (line). 



ties and to facilitate determination of minimum energy 
structures ("inherent structures"). Ref. [T^ provides 
a detailed specification of the modified potential. The 
modified potential is indistinguishable from the original 
BKS potential in terms of averaged structural and dy- 
namical behavior. As described in Ref. T5], the values of 
thermodynamic properties are slightly shifted compared 
to the original BKS potential, but the qualitative be- 
havior is unaffected. The Coulombic contribution to the 
energy is evaluated via the Ewald method, where the re- 
ciprocal space summation is carried out to a radius of 9 
times the smallest reciprocal cell width [l^- In all cases, 
the time step used in our MD simulations is 1 fs. 

We restrict our attention to the liquid phase, and the 
crystal phases stishovite, coesite and /3-quartz. A number 
of other crystal phases of silica are known. However, the 
stability fields of these three crystals dominate the phase 
diagram of silica on the widest scale of P and T, and are 
therefore natural first choices for examination 3] . These 
three crystals are also representative of the main types of 
local coordination structures found in silica crystals. The 
structure of /3-quartz is an open network of corner-shared 
Si04 tetrahedra; coesite is a denser network of corner- 
shared Si04 tetrahedra; and stishovite is a network of 
corner and edge shared SiOg octahedra. Previous work 
has shown that the BKS model is appropriate for study- 
ing both low and high density crystal structures [13 . To 
determine the phase diagram, our approach is to evalu- 
ate numerically the Gibbs free energy G of each of the 
phases as a function of P and T, and then seek the lines 
of intersection of these surface functions. 



A. Liquid free energy 

For the liquid phase, we use much of the equation of 
state data reported in Ref. [13, plus some new simula- 
tion data generated using the same methodology. These 
simulations modeled a system of a fixed number of 444 
molecular units (1332 ions) in the liquid phase along nine 
isochores from volumes V = 4.6296 to 8.6804 cm'^ mol~^, 
and ranging in T from nearly 7000 K to less than 2500 K 
(Fig. nj. Each of these liquid state points was equili- 
brated at constant V, and using velocity rescaling to at- 
tain a desired T. Average values of P and T were evalu- 
ated from subsequent constant- A'^Vi? runs having a du- 
ration of ten times the time required for silicon atoms to 
diffuse an average of 0.2 nm. The results provide E(T) 
and P(T) along the specified isochores. Ref. Il3| also 
describes the details of a calculation of the entropy of 
the liquid phase, Sr = 75.986 ±0.176 J mol-i K^^, at 
a reference state located at Tr = 4000 K and Vr = 
8.6804 cm3 mol'^ 

The value of E at an arbitrary point (Vo,To) on the 
surface E{V, T) is evaluated as follows. Along each of the 
nine isochores simulated, a 3rd order polynomial in T is 
fit to the E data. The value of E at the desired T = To 
is calculated from the polynomial found for each V . This 
creates a set of points approximating the curve E{V) at 
T — To- A cubic spline passing through these points is 
then found, creating a continuous function rcpresentating 
E{V) at T = To. The value of E at the point {Vo,To) 
is evaluated from this function. The value of P at an 
arbitrary point (Vo,To) is calculated in exactly the same 
way as for E, except that a 4th order polynomial in T 
is fit to the P data along each of the nine simulated 



isochores. Examples of the simulated and fitted values 
of E and P are shown in Fig. [3 

The value of S at an arbitrary point (Vo,To) is evalu- 
ated by thermodynamic integration, using the E{V,T) 
and P{V, T) surfaces constructed as described above. 
The integration is given by, 



-1 .8698 



S{Vo,To) 




P{V,To)dV 

(1) 

The definite integral over T is evaluated analytically us- 
ing the polynomial representation of i^ as a function of 
T, on the reference isochore. The definite integral over V 
is evaluated numerically via Simpson's rule, using data 
from the cubic spline representation of P as a function 
of V, constructed along the desired isotherm. 

We combine these numerical estimates to determine 
the Gibbs free energy G at arbitrary state points using 
G{V,T) = E{V,T) + VP{V,T) - TS{V,T). To find G 
at an arbitrary (P, T) point, we find the value of V from 
P{V, T) such that P has the desired value. In this way 
we can construct arbitrary isotherms or isobars cutting 
through the G{P, T) surface. 



B. Crystal free energy 

We conduct simulations of three crystal phases: 
stishovite, /3-quartz and coesite. Our simulations em- 
ploy 1200 ions for stishovite, 1536 ions for coesite, and 
1350 ions for /3-quartz. We carry out simulations over a 
range of V and T appropriate for each phase, as shown 
in Fig. n 

We employ the following procedure to obtain equilib- 
rium averages for E{T) and P{T) along the specified 
isochores. The rationale underlying this procedure is 
to allow us to obtain thermodynamic properties along 
a set of specified isochores, so that we may construct 
the G{P, T) surface for each phase in the same way as 
described above for the liquid phase. However, obtain- 
ing isochoric data for crystals requires care, as unit cell 
parameters may change with T, even though the over- 
all density remains fixed. In particular, we must ensure 
that anisotropic stresses do not arise in the simulation 
cell. The procedure, for each crystal, is as follows: 

First, we create an initial configuration of stishovite, 
/3-quartz {2(| and coesite 01. Then for a number of 
specified V , we optimize the T = atomic coordinates 
and unit cell parameters to minimize the energy and to 
remove anisotropic stresses. This optimization is car- 
ried out at constant overall V ^ and consists of alternat- 
ing applications of the simplex method (to optimize cell 
parameters) and the conjugate gradient method (to opti- 
mize atomic coordinates) |23| • This optimization cycle is 
repeated until the energy converges to a minimum value 
to within a tolerance of 10"^°. 

Then, for each T at which we desire thermodynamic 
properties, we carry out the following steps: 
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FIG. 3: Eanh{T) for (a) /?-quartz, (b) coesite and (c) 
stishovite along their respective reference isochores. In (a) the 
solid line is the fit to the data given by Eq.|21with Nmax = 2. 
The fit for Nmax = 3 is not visible as it overlaps with the 
Nmax = 2 curve on the scale of this plot. In (b) both the 
fit for Nmax ~ 2 (solid) and Nmax = 3 (dashed) are shown. 
In (c) the fits for Nmax ~ 2 (solid) and Nmax = 4 (dashed) 
are shown, while the curve for Nmax = 3 is not visible as it 
overlaps with the Nmax = 4 curve. 



(i) Beginning with the optimized configuration at the 
appropriate V, we conduct a 20 ps constant V simulation, 
during which the desired T is established via velocity 
rescaling every 100 time steps. 

(ii) The configuration produced in (i) is used to initi- 
ate a 20 ps constant-iVy_E simulation, to ensure that an 
equilibrium state at the desired T has been achieved. 

(iii) To relax any anisotropic stress that may have 
arisen in bringing the system to non-zero T, we carry 
out a 40 ps constant- A^PT simulation (during which the 
simulation cell geometry is unconstrained) where we set 
P to the average value from step (ii). 

(iv) Step (iii) may have changed the overall V of the 
simulation cell away from the desired isochore. We re- 
store the value of V of the simulation cell by isotropi- 
cally rescaling the average cell lengths obtained in step 
(iii), while leaving the obtained average angles fixed. For 
all crystals, the rescaling is never more than 0.5% of the 



desired volume, and is typically 0.1%. This rescaled con- 
figuration is then used to initiate a 30 ps constant-iVT^i? 
simulation, during which the average values of P and T 
are evaluated. 

We note that for /3-quartz, step (iii) is carried out for 
50 ps and step (iv) for 80 ps. These longer times are used 
in order to resolve the subtle variation of P along iso- 
chores, since /?-quartz displays a density maximum in the 
region of our simulations. We also note that the /3-quartz 
phase spontaneously converts to a-quartz, but only for 
T and V outside the range of simulated points shown for 
/3-quartz in Fig. ^ Our results therefore pertain only to 
/3-quartz and are not influenced by this crystal-crystal 
phase transition. 

The above procedure provides E{T) and P{T) along 
specified isochores. Using the same fitting and interpola- 
tion procedure as is used for the liquid, we can therefore 
evaluate E and P at arbitrary state points {V,T). 

Finally, we need to evaluate S for each crystal at a 
reference state point, in order to construct the surface 
S{V,T) via thermodynamic integration. Our method is 
as follows. For each crystal phase we select a reference 
volume Vu (see Fig.^, and choose the Tr = 1500 K con- 
figuration obtained at the end of step (iv) above. Using 
the conjugate gradient method, we optimize the atomic 
positions (at fixed cell geometry) to find the minimum 
energy configuration. We then evaluate the Hessian ma- 
trix of this minimum energy configuration and diagonal- 
ize it to find the eigenfrequency spectrum. The classical 
harmonic entropy is found from this eigenfrequency spec- 
trum. (The details of this approach are given in Ref. [Tjl , 
where the method is used to find the classical harmonic 
entropy of inherent structures of the liquid state.) 

To determine the total entropy, we need to evaluate 
the anharmonic contribution and add it to the harmonic 
entropy found above. We use the energy-optimized con- 
figuration for which we calculate the harmonic entropy 
as the starting configuration for 15 equally spaced sim- 
ulations from T = 100 K to Tr = 1500 K. We simulate 
each state point at constant V using velocity scaling to 
maintain T at the desired value, with fixed cell geometry, 
for 150 ps (400 ps for stishovite). From these simulations 
we evaluate. 



EanhiT) = U{T) - -R{\ - 1/N)T, 



(2) 



where U is the potential energy and R is the gas constant. 
Using a polynomial fit, 
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we evaluate. 
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FIG. 4: AG, the Gibbs free energy relative to tliat of coesite 
(C) for the hquid (L), /3-quartz (Q) and stishovite (S) phases, 
at constant P — 2 GPa. The intersections locate points on the 
stable and metastable coexistence lines that cross this isobar. 



we obtain error estimates for Sanh{Tfi)- Fig. O 
shows the variation of Eanh with T for each of the 
three crystals simulated. The resulting reference en- 
tropies for each crystal phase at Tr — 1500 K at 
their respective reference volumes Vr are as follows: 
44.982 J mol-i K'^ at Vr = 8.9933 cm^ mol'^ for /3- 
quartz; 43.682 J mol'^ K'^ at Vr = 7.1478 cm^ mor^ 
for coesite; and 39.536 J mol~^ K~^ at Vr — 
4.7650 cm"^ mol^^ for stishovite. The uncertainty in each 
S value is approximately 0.01 J mol""'^ K"-*^. 

The above procedure provides E{V,T), P{V,T) and 
the reference value of S for each of the three crystal 
phases. The procedure to evaluate G{P, T) from this 
information for each crystal phase is the same as is used 
for the liquid phase. 



C. Coexistence boundaries and error estimates 

For every pair of phases we determine the coexistence 
line as the locus of points in the plane of P and T for 
which G for the two phases is the same. Along each locus, 
we also find the value of V for each of the two coexisting 
phases. Fig. 0] shows an example of the intersection of 
isobars of G (relative to coesite) for each phase at P = 
2 GPa. 

We perform several checks on our scheme. We calcu- 
late the change in S for a single phase around a closed 
path in the V-T plane, which we find to be zero within 
an error of approximately 0.01 J mol^^ K~^. We also 
check that the relations P — ~{dA/dV)T (where A 
is the Helmhohz free energy) and P ~ T{dP/dT)v = 
— {dE/dV)T are satisfactorily met. Furthermore, along 
the coexistence lines, we check the Clapeyron relation 
dP/dT = AS/AV, where AS is the difference in S be- 
tween the two phases and AV is the difference in V; we 
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FIG. 5: (a) Experimentally determined coexistence lines of 
silica in the P-T plane. Stability fields for the stishovite (S), 
coesite (C), /3-quartz (Q) and liquid (L) phases are shown. 
Both stable (solid) and metastable (dashed) coexistence lines 
are shown. The inset shows the stability fields of cristobalite 
and tridymite, not considered in this work. Adapted from 
Ref. Q. (b) Phase diagram of BKS silica in the P-T plane. 
Solid lines are stable coexistence lines. Dotted lines show error 
estimates for the crystal-liquid coexistence lines, as described 
in the text. Metastable coexistence lines (dashed) are also 
shown that meet at the metastable S-L-Q triple point. 



find this to be satisfied to within 0.1 MPa K~^. 

Throughout the evaluation scheme described above, 
the largest single source of statistical error is the uncer- 
tainty cited in Ref. jl^ for Sr, the entropy of the liquid 
at the reference state point. We therefore create confi- 
dence limits for our melting hnes, shown in Fig. [5] by 
allowing the value of Sr to vary by ±0.18 J mol^^ K^^. 



III. RESULTS AND DISCUSSION 

Fig. intb) plots P-T coexistence conditions, both sta- 
ble and metastable, occurring among the liquid phase 




V (cm mol ) 



FIG. 6: Phase diagram of BKS silica in the V-T plane. The 
notation and symbols used have the same meaning as in Fig.|S] 
Note that in this projection, both one-phase stability fields as 
well as two-phase coexistence regions are located. The pro- 
jections of the metastable coexistence lines (dashed) shown in 
Fig. 13 are also presented. 



(L) and the crystalline phases /3-quartz (Q), coesite (C) 
and stishovite (S). Fig. |H| is the projection of the same 
boundaries onto the plane of V and T. This plot ex- 
poses the volume differences of coexisting phases along 
phase boundaries. This type of plot is rarely constructed 
for real materials, due to the challenge of determining the 
densities of coexisting phases, especially at high pressure. 
However, it is readily constructed from simulation data. 

Comparison of the BKS and experimental phase 
boundaries Ql in Fig.|31exposes the quantitative deficien- 
cies of the model. Apparent in particular is the difference 
between the pressures at which corresponding features 
occur. For example, the S-L-C triple point occurs at 
13.4 GPa in real silica, but at only 5.8 GPa in the model. 
Overall, the P range of the crystal stability fields is sub- 
stantially lower in the model. The correspondence of the 
thermal behavior is better, but significant differences still 
occur. The T of the S-L-C and C-L-Q triple points are 
respectively 15% and 32% higher than their experimental 
values. 

However, it is noteworthy that the topology of the real 
silica phase diagram is reproduced by the model. All 
three studied crystals have large stability fields, which 
increase in extent as P increases. More subtle features, 
notably the melting line maxima in both the Q-L and C- 
L coexistence lines are also reproduced. Also, the occur- 
rence of a metastable S-Q-L triple point in the stability 
field of coesite, suggested by an extrapolation of the ex- 
perimental boundaries, is observed in the model. Thus, 
despite its quantitative deficiencies, the BKS model is 
appropriate for studying the qualitative behavior of a 
substance with a silica-like phase diagram. 

The phase information given in Figs. [S] and |3 allows 
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previous (and future) observations of the behavior of 
BKS silica to be considered within the context of the 
phase behavior of the model itself. For example, sev- 
eral studies of BKS silica have identified the location of 
a density maximum in the liquid phase ^, 13] . A ther- 
mal anomaly, a line of maxima of the isochoric specific 
heat Cy, has also been located in simulations of the liq- 
uid [ij, [13 . Both of these features have been related to 
the early stages of the formation of a structured tetra- 
hedral network in the liquid state. This structural evo- 
lution also is believed to underlie a crossover from non- 
Arrhcnian ("fragile") to Arrhenian ("strong") dynamics 
in the liquid ,13l.14.15j . 

We show in Fig. [7| the location of the line of density 
and Cy maxima in both the P-T and V-T planes. These 
lines approximately separate the liquid behavior into a 
"tetrahedral network influenced" region at low T and 
large V (low P), and a "normal liquid" region at high 
T and small V (high P). Consistent with this, the sta- 
bility fields of /3-quartz and coesite (both of which have 
four-coordinated silicon atoms) occur within the network 
influenced region, while the stability field of stishovite 
(with six-coordinated silicons) falls outside. 

We also show in Fig. [T] the location of the state point 
at which evidence of liquid-liquid phase separation was 
reported in Ref. 5]. This point occurs at a density just 
above that of the high density edge of the one-phase sta- 
bility field of coesite. This is a plausible density at which 
the open tetrahedral network structure of the (one-phase) 
supercooled liquid state begins to collapse to a higher 
density, perhaps via a discontinuous phase transition. 



FIG. 7: BKS phase boundaries in (a) the P-T plane and (b) 
the V-T plane, in relation to density maxima (filled circles) 
and Cv maxima (squares) in the liquid phase. Also located is 
the state point (star) at which evidence of liquid-liquid phase 
separation was reported in Ref. Q- 
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